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Abstract —In this paper, least square estimation (LSE)-based 
dynamic generator model parameter identification is investigated. 
Electromechanical dynamics related parameters such as inertia 
constant and primary frequency control droop for a synchronous 
generator are estimated using Phasor Measurement Unit (PMU) 
data obtained at the generator terminal bus. The key idea of 
applying LSE for dynamic parameter estimation is to have a 
discrete autoregression with exogenous input (ARX) model. With 
an ARX model, a linear estimation problem can be formulated 
and the parameters of the ARX model can be found. This paper 
gives the detailed derivation of converting a generator model with 
primary frequency control into an ARX model. The generator 
parameters will be recovered from the estimated ARX model 
parameters afterwards. Two types of conversion methods are 
presented: zero-order hold (ZOH) method and Tustin method. 
Numerical results are presented to illustrate the proposed LSE 
application in dynamic system parameter identification using 
PMU data. 

Index Terms —Least squares estimation (LSE), Phasor Mea¬ 
surement Unit (PMU) 

1. INTRODUCTION 

Traditionally, Supervisory Control and Data Acquisition 
(SCADA) system using nonsynchronous data with low density 
sampling rate is used for monitoring and control of the system. 
Such measurements can not capture the system dynamics. The 
advent of phasor measurement units (PMUs) equipped with 
GPS antenna provides voltage/current phasors and frequency 
with a high density sampling rate up to 60 Hz. These phasor 
measurements transmitted with time stamps can help control 
systems have an accurate picture of the power system. 

Synchronous generator parameter estimation has been in¬ 
vestigated in the literature. Based on the scope of estimation, 
some only investigate electrical state estimation (e.g. rotor 
angle and rotor speed) ||T|, Q, while others estimate both 
system states and generator parameters 0 0 Based on 
estimation methods, there are at least two major systematic 
methods for parameter estimation: least squares estimation 
(LSE) ||7j-|9|| and Kalman filter-based estimation (D-GD- 
To use LSE for dynamic system parameter estimation, a 
window of data is required. On the other hand, Kalman filter- 
based estimation can carry out estimation procedures at each 
time step. Thus Kalman filter-based estimation can be used 
for online estimation. 

Majority of the research papers on PMU-based dynamic 
parameter estimation adopt Kalman filter-based estimation 
approach, e.g., |T0|- OD- In the literature, there are plenty of 
research on LSE-based offline generator parameter estimation 
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based on either time-domain data or frequency response data 
0-0 GD -(T^. However, little research can be found 
for LSE-based dynamic parameter estimation for PMU data- 
based applications. In addition, existing research on LSE- 
based generator parameter estimation uses nonlinear LSE 0. 
existing toolboxes for nonlinear system identification (e.g., 
Matlab System Identification Toolbox) The adoption of 
Matlab toolbox does not provide the insights as how the 
algorithms work, while the nonlinear LSE often time faces 
convergence issues if the set of the parameters chosen are not 
suitably chosen as demonstrated in Q. 

Therefore, this paper aims to adopt linear LSE for generator 
parameter estimation via PMU data. Electromechanical dy¬ 
namics related parameters such as inertia constant and primary 
frequency control droop for a synchronous generator are 
estimated using PMU data obtained at the generator terminal 
bus. The key idea of applying LSE for dynamic parameter 
estimation is to have a discrete ARX model. With an ARX 
model, a linear estimation problem can be formulated and the 
parameters of the ARX model can be found. This paper gives 
the detailed derivation of converting a generator model with 
primary frequency control into an ARX model. The generator 
parameters will be recovered from the estimated ARX model 
parameters afterwards. Two types of conversion methods are 
presented: zero-order hold (ZOH) method and Tustin method. 
Numerical results are presented to illustrate the proposed LSE 
application in dynamic system parameter identification using 
PMU data. 

The rest of the paper is organized as follows. Section II 
describes the LSE for an ARX model. Section III describes 
the two approaches to convert a continuous transfer function 
into a discrete ARX model. Section IV presents the step-by- 
step procedures from a continuous generator model to ARX 
models along with numerical estimation results. Section VI 
concludes the paper. 

H. LSE FOR ARX Model 

A block diagram of a simplified synchronous generator 
model is shown in Eigure where S is the rotor angle (rad), 
cjQ is the synchronous speed (rad/s), U is the inertia constant 
(seconds), D is the damping factor, R is the speed regulation 
constant (p.u.), T is the time constant of the turbine-governor 
(seconds), Pm is the mechanical power input (p.u.), and Pe is 
the electrical power (p.u.). 

Linearizing the above model, we can develop the 
Continuous-Time (CT) closed-loop transfer function in 
Laplace domain, which will be converted into its equivalent 
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Fig. 1. A simplified synchronous generator model. 


Discrete-Time (DT) transfer function in Z-domain where sam¬ 
pling is considered. The z-transfer function is then converted 
into an autoregressive exogenous (ARX) model as follows: 
Suppose that the DT transfer function is given by: 

yjz) _ biz~'^ + boz-"^ 

u{z) z ^ aiz~^ ^ aoz~‘^ 

where y is the output, and u is the input. Then, its equivalent 
ARX model can be expressed as: 

y{k) = —aiy{k—l)—af^y{k—2)-^hiu{k—l)^h{)u{k—2)^e{k) 

( 2 ) 

where e is the error. 

Given time series measurements of the input and o^ut, 
an overestimated problem can be formulated based on S as 
follows: 
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In a concise format, can be written as 6 = Ax + e. The 
parameter vector x can be found as {AA^)~^b. 

Converting the generator’s dynamic model into a difference 
equation, LSE-based estimation can then be used to estimate 
the coefficients. After that, the parameters H, D, T, and R can 
be identified. Before proceeding any further, let’s consider how 
to discretize a CT transfer function given in Laplace domain 
to find its equivalent 2 ;-transfer function. 

III. DISCRETIZATION OF A CT TRANSFER FUNCTION 

Several methods exist for discretizing a CT transfer function 
given in Laplace Domain. These methods include: 1) Zero- 
order hold (ZOH) method, where a zero-order hold element 
is placed at the input of the system to hold the input signal 
constant during each sampling interval. 2) Numerical approxi¬ 
mations to time integrals, where Euler’s forward, Euler’s back¬ 
ward, or Tustin’s approximations are used to get substitution 
formulas for the Laplace operator s in terms of z and the 
sampling period h. Although the latter methods are much 
easier to compute, care should be taken with Euler’s forward/ 
backward approximations because they cannot be used without 
modifications as discussed in p0| . Eor this reason, only ZOH 
and Tustin’s approximation methods are investigated in this 
paper. 

A. Zero-order hold method 

Eor a general linear system, the relationship between the 
input and output can be determined from the response of the 
system to a given signal, such as a unit-step input pT| . Let 


H{s) denote the CT transfer function from the input u{t) to 
the output y{t). Assuming that the system is preceded by a 
zero-order hold element and followed by a sampler as shown 
in Eigure we wish to find the DT transfer function from the 
sampled input u{k) to the sampled output y{k). This can be 
achieved by determining the step response of the CT system 
first. After that, the corresponding 2 ;-transform of the step 
response is obtained. Einally, we divide by the 2 ;-transform 
of the unit-step function. 



H{z) 


Fig. 2. Discretization of a continuous-time system. 


This procedure can be formulated as follows: 



where is the ^-transform, ^ ^ is the inverse Laplace, and 
h is the sampling interval. 

R.^Tustin's method 

Numerical approximations to the time integral (or deriva¬ 
tive) are widely used in digital control applications. [REE 
2books]. A DT transfer function in z-domain is obtained by 
using appropriate substitutions to time integrals so that the 
behaviors of the z-transfer function and the CT s-transfer 
function are similar. One of these methods is Tustin’s approx¬ 
imation which is also known as the bilinear transformation. 

In order to obtain the substitution formula for Tustin’s 
method, consider the simple integrator y = u, where u is 
the input and y is the output. In Laplace domain the above 
equation can be written as: = 2. 

On the other hand, Tustin’s method uses the following 
approximation to the time integral: 

y{k) « y(/c - 1) + ^ [u{k) + u{k - 1)] (5) 

In z-domain, the difference equation 0 corresponds to: 
y{z) ^,h z + l 

u{z) ^ 2 z-1 ^ ’ 

We finally find the substitution formula: s •(— f 

As a result, Tustin’s approximation maps the left-half of the 
s-domain into the unit circle |z| < 1 in z-domain as shown in 
Eigure 


IV. Case studies 

Eor the case studies, a benchmark model based on Eigure 
is built using MATLAB/Simulink to generate the data where a 
step response in Pg is simulated. The input and output data are 
collected and then fed into the estimation algorithm to find out 
the coefficients and bi. After that, the parameters P, P, T, 
and D are recovered for several sampling intervals h. Different 












































3 



Fig. 3. Mapping of the stability region for Tustin’s approximation. 

scenarios are simulated and the results are presented for each 
case. The following parameters are used in the simulation: 
H = 2.5, R = 0.05, and T = 0.5. 


If we can estimate the coefficients ai, uq, and ho, we 

can find out the parameters. The ARX model can be put 

in the matrix format b = Ax as follows: 
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Using LSE, the above overestimated problem can be solved 
hy X = A)~^A^h. 

After solving the optimization problem and finding the 
coefficients, the parameters can be found using the following 
equations: 


A. Case 1: Power (APe) versus frequency (Auj) 

In order to illustrate the procedure and to simplify the 
computations, the damping factor is ignored (i.e., D = 0) 
in this case. 

1) ZOH method: Consider the synchronous generator’s 
model shown in Figure Without the damping factor, the 
CT closed-loop transfer function from APe to Auj is given 
by: 


H{s) = 


Auj 

'aK 


TsPl 


2HTs^ 


(7) 


Since ^ > H, the transfer function represents an under¬ 
damped second-order system with two complex conjugate 
poles given by: 


Sl,2 = ^±j 




2HT 


( 8 ) 


The transfer function H{s)/s will be fractionated and then 
inverse Laplace transformation will be applied. 

It worth mentioning that even if we include the damping 
factor, say D = 0.8 p.u., we still get an underdamped 
response which is generally the case for the synchronous 
generator under normal operating conditions as discussed in 
Following the procedure given by 0 , we obtain the DT 
transfer function as: 

y{z) hz + bo 


H{z) 


(9) 


u{z) z^ + aiz + ao 

Thus, the ARX model can be obtained based on 0 as follows: 

y(k-h2) = -aiy(k-hl)-aoy(k)-hbiu(k-hl)-hbou(k) (10) 
where y is the output Auj and u is the input APg, 

= — 2 e^ cos {ujh) 


ai 

ao 

^0 
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= e T 


= P ( 1 — e 2 t cos (ujh) — ke sin 
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( 11 ) 


e T — e 2T cos {ujh) 


ke^T sin 


and a, uj, and k are given by the following equations: 
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2) Tustin's method: Obtaining the DT transfer function 
using Tustin’s approximation method is much easier than the 
ZOH method because we do not need to worry about the 
poles or the system’s type of response. Substitute the Laplace 
operator s in we get: 

= (15) 

u[z) + aiz + ao 


where the coefficients are given by: 
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ao 
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2-AHRTk‘^ 
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= ^ 
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(16) 


and a = 2HRTkf + 2HRk + 1, /c = |. After that, the 
ARX model is developed and the optimization problem is 
solved as discussed before. Once the coefficients are found, 
the parameters can be easily found. Expressions of bi and 62 
are used to find T as follows: 


T 

B 


_ 2B-1 
k ’ 



(17) 


then, using the expression of ai we can find R in terms of H 
as: 


R = 


2 — ai 

i7(2aiTP + 2aiA: + 4TP) 


(18) 


Finally, we can substitute •dll' and ( [Tsl l in the expression of 
62 to find H as follows: 


1 + TA: - ^ 

2b2k{Tk + 1 ) 


(19) 
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3) Numerical results: The benchmark model is simulated 
using MATLAB/Simulink for various input signals and the 
parameters are estimated. Noise with 0.0001 variance and zero 
mean is added to the input signal. A step input equals 0.2 p.u. 
is applied at t = 1 second and the p.u. change in the input 
power (APe) and output frequency (Acc) are shown in Figure 
IHwhen h = 0.1 second. 



_0.02-i-i-i-1-1-i-1-1-i- 

0123456789 10 
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Fig. 5. Change in power APe and corresponding change in angle Ad. 


1) ZOH method: Following the aforementioned procedure, 
the DT transfer function has the form: 

h2Z^ + biz + 6 o 


H{z) = 


yjz) ^ _ 

u{z) + a 2 z'^ + aiz + uq 


( 20 ) 


Fig. 4. Change in power APe and corresponding change in frequency Auj. 


The estimated coefficients using Matlab Optimization toolbox 
CVX when h = 0.1 are shown in Table |l| And the re¬ 
covered parameters for different sampling intervals are shown 
in Table ini 

TABLE I 

ESTIMATED COEEEICIENTS USING CVX 


ine Coefficient 
ine 62 
ine bi 
ine bo 
ine a I 
ine ao 
ine 


ZOH method 

X 

19.7333x10-3 

-16.1380x10-3 

-1.7467 

818.5742x10-3 


Tustin’s method 
9.8214x10-3 
1.7857x10-3 
-8.0357x10-3 
-1.7500 

821.4286x10-3 


TABLE II 

Recovered parameters eor dieeerent h 


ine 


ZOH method 

Tustin’s method 

ine 

h= 0.1 

h= 0.01 

h= 0.001 

h= 0.1 

h= 0.01 

h= 0.001 

ine T 

0.499 

0.499 

0.500 

0.500 

0.500 

0.500 

ine R 

0.050 

0.050 

0.050 

0.0500 

0.0499 

0.0499 

ine H 
ine 

2.506 

2.500 

2.499 

2.500 

2.499 

2.499 


It should be mentioned that in order for Tustin’s method to 
work, the CT model in Laplace domain cannot be simulated 
directly, otherwise, the estimated coefficients are almost equal 
to the coefficients of ZOH method with 62 = 0- Discrete 
transfer function blocks shall be used where each s is replaced 
by this way correct values of the parameters are 

obtained. 


where: 


' 02 


ai 

ao 

b2 

hi 

ho 


= —2e 2 t cos {(jjh) — 1 
= + 2e^ cos {ujh) 

-h 

= —e T 


= R 
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( 21 ) 


and a = T-2HR, b= , and uj = ^• 

After estimating the coefficients, the parameters can be 
identified. The recovered parameters for different sampling 
intervals are listed in Table [nil 

2 ) Tustin ’s method: The DT transfer function is in the form: 


where: 


TJ! ^ y{z) bzz^ + b2z‘^ + biz + bo 
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u[z) z'^ + a 2 Z^ + aiz + ao 


02 

Ol 
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< 63 
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^ -2HRTk^+2HRk‘^-k 
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( 22 ) 


(23) 


and k = a = 2HRTk^ + 2HRk^ + k. The recovered 
parameters after estimating the coefficients are shown in Table 

mil 


B. Case 2: Power (AP^) versus rotor angle (AS) 

In this case, the change in the rotor angle A^ is considered 
as the output while the change in the electrical power APg is 
the input. The damping factor is again ignored in this case. 
The same input of case 1 is applied at t = 1 second and the 
p.u. change in both power and angle are shown in Figure 
for h = 0.1 second. Since the procedure is explained in the 
first case, it will not be covered here. 


TABLE III 

Recovered parameters using ZOH and Tustin’s methods 


ine 


ZOH method 

Tustin’s method 

ine 

h= 0.1 

h= 0.01 

h= 0.001 

h= 0.1 

h= 0.01 

h= 0.001 

ine T 

0.499 

0.499 

0.499 

0.499 

0.500 

0.500 

ine R 

0.050 

0.050 

0.0499 

0.050 

0.0499 

0.0498 

ine H 
ine 

2.501 

2.500 

2.499 

2.499 

2.500 

2.500 
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C. Case 3: Real-world PMU data 

In this case both ZOH and Tustin’s methods, are used to 
estimate the parameters from a recorded real-world PMU data 
of an anonymous busbar of MISO system. The data is recorded 
for about 40 seconds when a major disturbance occurs in the 
system. Similar to case 1, the change in power and frequency 
are used for estimating the parameters. The estimated param¬ 
eters for both methods are listed in Table |IV| which clearly 
shows that Tustin’s method fails to provide correct estimation 
of the parameters. In order to validate the results of ZOH 
method, the model is built in MATLAB/Simulink and event 
playback is used to inject the change in power as input. The 
change in frequncy is then compared to the PMU measurement 
as shown in Figure which demonstrates a great degree of 
match. 


TABLE IV 

Recovered parameters eor real-world PMU data 


ine 

T 

R 

H 

ine ZOH method 

0.4534 

0.2320 

13.8945 

ine Tustin’s method 

0.001 

0.000667 

0.0639 

ine 






Time (sec) 


Fig. 6. Real-world PMU measurements versus ZOH estimation. 


V. CONCLUSION 

In this paper, linear LSE-based generator model parameter 
estimation via PMU data has been presented in detail. The 
continuous time Laplace model is first converted to a discrete 
ARX model using ZOH method or Tustin’s method. The 
coefficients of the ARX model will then be identified using 
linear LSE. The generator model parameters are then recov¬ 
ered by the ARX model coefficients. The proposed approaches 
are illustrated in numerical case studies to demonstrate their 
effectiveness in identifying generator parameters. It is found 
that for real-world data, ZOH method is robust while Tustin’s 
method is very sensitive to the estimation model. 
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